

%% load everything of interest
load('data_direction.mat')
load('sim_eq_direction_round1.mat')
load('parameter_estimates_round1_replication.mat','specif')
[DMES,intervEq] = createDMES(intervEq,V_R,V_US,V_UK,V_France,V_Russia,V_China,YY,specif);


%% normalize continuous variables
% data.gdp_pc = (data.gdp_pc-mean(data.gdp_pc))/std(data.gdp_pc);
% data.pop = (data.pop-mean(data.pop))/std(data.pop);
% data.polity = (data.polity-mean(data.polity))/std(data.polity);
data.terrain = (data.terrain-mean(data.terrain))/std(data.terrain);
dd = [data.USdist,data.UKdist,data.FRdist,data.RUdist,data.CHdist];
dd = (dd-mean(dd(:)))/std(dd(:));
data.USdist = dd(:,1);
data.UKdist = dd(:,2);
data.FRdist = dd(:,3);
data.RUdist = dd(:,4);
data.CHdist = dd(:,5);
clear dd


%% estimates of mean utilities
SV=25; % # of starting values
X=[ones(N,1) data.terrain]; K=size(X,2);
Z_US=[data.USdist data.USally data.UScolony data.USwar]; KZ=size(Z_US,2);
Z_UK=[data.UKdist data.UKally data.UKcolony data.UKwar];
Z_France=[data.FRdist data.FRally data.FRcolony data.FRwar];
Z_Russia=[data.RUdist data.RUally data.RUcolony data.RUwar];
Z_China=[data.CHdist data.CHally data.CHcolony data.CHwar];
Z=[Z_US,Z_UK,Z_France,Z_Russia,Z_China];
reb={1:K,2:KZ,2:KZ};
mpow={[1,3:K],[1,3:K],1,1};
Y=arrayfun(@(n)find(sum(YY==table2array(outcomes(n,4:end)),2)==size(YY,2)),1:N)';

ltheta=length(reb{1})+4+length(mpow{1})+4+length(mpow{2})+5+length(reb{2})+5+...
    length(reb{3})+length(mpow{3})+length(mpow{4})+20+length(specif);

% coefficient restrictions
A=[[zeros(4,length(reb{1})+2*(4+length(mpow{1}))),...
    [1,-1,0,0,0;1,0,-1,0,0;1,0,0,-1,0;1,0,0,0,-1],...
    zeros(4,ltheta-length(reb{1})-2*(4+length(mpow{1}))-5)]; % gamma_m^1=gamma_{m'}^1
    [zeros(4,length(reb{1})+2*(4+length(mpow{1}))+5+KZ-1),...
    [1,-1,0,0,0;1,0,-1,0,0;1,0,0,-1,0;1,0,0,0,-1],...
    zeros(4,ltheta-length(reb{1})-2*(4+length(mpow{1}))-2*5-(KZ-1))]; % gamma_m^2=gamma_{m'}^2
%     [zeros(KZ-1,length(reb{1})+2*(4+length(mpow{1}))+5),eye(KZ-1),...
%     zeros(KZ-1,length(thetaStarExp)-length(reb{1})-2*(4+length(mpow{1}))-5-(KZ-1))]; % gamma_0^1=0
%     [zeros(KZ-1,length(reb{1})+2*(4+length(mpow{1}))+5+KZ-1+5),eye(KZ-1),...
%     zeros(KZ-1,length(thetaStarExp)-length(reb{1})-2*(4+length(mpow{1}))-2*(5+(KZ-1)))]; % gamma_0^2=0
    [zeros(4+length(mpow{1}),length(reb{1})),eye(4+length(mpow{1})),-eye(4+length(mpow{1})),...
    zeros(4+length(mpow{1}),ltheta-length(reb{1})-2*(4+length(mpow{1})))]; % phi_m^1=phi_m^2
    [zeros(14,ltheta-length(specif)-20),...
    [1,-1,zeros(1,8),zeros(1,10);
    1,zeros(1,3),-1,zeros(1,5),zeros(1,10);
    zeros(1,2),1,-1,zeros(1,6),zeros(1,10);
    zeros(1,2),1,zeros(1,2),-1,zeros(1,4),zeros(1,10);
    zeros(1,2),1,zeros(1,3),-1,zeros(1,3),zeros(1,10);
    zeros(1,2),1,zeros(1,4),-1,zeros(1,2),zeros(1,10);
    zeros(1,2),1,zeros(1,5),-1,0,zeros(1,10);
    zeros(1,10),1,-1,zeros(1,8);
    zeros(1,10),1,zeros(1,3),-1,zeros(1,5);
    zeros(1,10),zeros(1,2),1,-1,zeros(1,6);
    zeros(1,10),zeros(1,2),1,zeros(1,2),-1,zeros(1,4);
    zeros(1,10),zeros(1,2),1,zeros(1,3),-1,zeros(1,3);
    zeros(1,10),zeros(1,2),1,zeros(1,4),-1,zeros(1,2);
    zeros(1,10),zeros(1,2),1,zeros(1,5),-1,0],...
    zeros(14,length(specif))]]; % "West" and "East" deltas

options=optimset('Display','off','HessFcn',...
    @(x,lambda)HessPhat(x,lambda,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow));

load('parameter_estimates_direction_round1_replication.mat','srngSV')
rng(srngSV)
theta0=randn(ltheta,SV)/10; % starting values
theta=theta0; 
l0=1:SV; % maximized likelihood values across starting values
exitn=1:SV; % Knitro exit flags across starting values
f0 = exitn;

time=0;
for sv=1:SV
    tic
    f0(sv) = Phat(theta0(:,sv),Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow);
    if ~isfinite(f0(sv))
        theta(:,sv) = NaN(length(theta0(:,sv)),1);
        exitn(sv) = NaN;
        l0(sv) = NaN; 
        disp(sv)
        toc
        time=time+toc;
    else
        [theta(:,sv),l0(sv),exitn(sv)]=knitromatlab(@(x)Phat(x,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow),...
            theta0(:,sv),[],[],A,zeros(size(A,1),1),[],[],[],[],options,'knitro1.opt');
        disp(sv)
        toc
        time=time+toc;
    end
end


%%
clearvars -except exitn K l0 N SV theta theta0 time X Y YY Z reb mpow srngSV specif
save('parameter_estimates_direction_round1.mat')


%% verify optimality
[~,use] = min(l0);
thetaStar = theta(:,use);
load('parameter_estimates_round1_replication.mat')
[~,use] = min(l0);
thetaStar0 = theta(:,use);
disp(' ')
disp('verify round 1 direction results:')
disp(['max. rel. diff. between estimates = ',num2str(max(abs(thetaStar0-thetaStar)./abs(thetaStar0)))])



